This week, we’ll be using a Java application called OpenTripPlanner to create isochrones (walksheds and drivesheds) from within R, using the opentripplanner package.

I’ll be using these libraries:

library(osmdata)
library(opentripplanner)
library(tidyverse)
library(sf)
library(ggthemes)
library(ggspatial)

You’ll also need to have Java installed on your computer, since we’ll be using a Java application called Open Trip Planner.

Load locations

I’m going to start by importing a KML file from the City of Cambridge’s Open Data Portal. This data set shows the locations of the all Cambridge Public Library locations. You can see from the results messages that this is point data with latitude/longitude coordinates, and we’ll keep them that way.

CPL_libraries <- st_read(
  "https://data.cambridgema.gov/api/geospatial/kn2z-b6eg?method=export&format=KML")
## Reading layer `Layer #0' from data source `https://data.cambridgema.gov/api/geospatial/kn2z-b6eg?method=export&format=KML' using driver `KML'
## Simple feature collection with 7 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -71.14684 ymin: 42.36391 xmax: -71.08483 ymax: 42.39266
## geographic CRS: WGS 84

Get street data

We’ll use Open Trip Planner to find the areas of Cambridge that are within five minutes of a library by walking, by car, and by bike. To start, we need data on the street network. Before doing that, you’ll need to do some set-up.

Create a folder in your Repo called OTP. Within that folder, create another folder called graphs. Within that folder, create another folder called default. Add the OTP directory to your .gitignore file.

Next, you can run the following code (but for the area you’re analyzing) to download street network data from Open Street Map.You only need to do this once. Once you have the *.osm folder in your OTP/graphs/default directory, you can delete this code chunk.

opq(bbox = 'Cambridge MA USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_xml(file = 'OTP/graphs/default/cambridge_streets.osm')

Optionally, you can also get a set of sf features from OpenStreetMap to plot on a map.

MA_state_plane <- "+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +ellps=GRS80 +units=m +no_defs"

cambridge_street_features <- opq(bbox = 'Cambridge MA USA') %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_sf()

cambridge_streets <- cambridge_street_features$osm_lines %>%
  st_transform(crs = MA_state_plane)

Verify that you now have a file in your default directory called “cambridge_streets.osm”, and try plotting the streets you’ve downloaded. You’ll notice that it isn’t clipped to the town boundaries. It includes all the roads in the rectangle (bounding box) that contains the town boundaries.

ggplot(cambridge_streets) +
  geom_sf() +
  theme_map()

Set up Open Trip Planner

Now we need to download a little Java utility called otp.jar and save it to the OPT directory we’ve created. You only need to run this line once, so you can delete this chunk after you run in - once you confirm that you have a file called otp.jar in your OPT directory on your computer

path_otp <- otp_dl_jar("OTP")
## The OTP will be saved to OTP/otp.jar

Now, we’ll build a graph. This is a representation of the street and transit networks.

path_data <- file.path(getwd(), "OTP")
path_otp <- paste(path_data, "otp.jar",sep = "/")

otp_build_graph(otp = path_otp, dir = path_data, memory = 1024)

And now we’ll launch Open Trip Planner! otp_setup may take several minutes to run.

After you do this next part, an OTP application will open in you web browser. This indicates that OpenTripPlanner is running, but other than that, you can ignore it.

otp_setup(otp = path_otp, dir = path_data, memory =1024)
## 2020-09-30 11:57:00 OTP is loading and may take a while to be useable
## Router http://localhost:8080/otp/routers/default exists
## 2020-09-30 11:58:00 OTP is ready to use Go to localhost:8080 in your browser to view the OTP
# Connect to opentripplanner
otpcon <- otp_connect()
## Router http://localhost:8080/otp/routers/default exists

Create isochrones

Create isochrones for areas within a five-minute walk and a five-minute drive.

iso_5min_walk <- 
  otp_isochrone(otpcon = otpcon, fromPlace = CPL_libraries, 
                mode = "WALK", cutoffSec = 300) %>%
  st_transform(crs = MA_state_plane) %>%
  mutate(mode = "walk")

iso_5min_drive <- 
  otp_isochrone(otpcon = otpcon, fromPlace = CPL_libraries, 
                mode = "CAR", cutoffSec = 300) %>%
  st_transform(crs = MA_state_plane) %>%
  mutate(mode = "drive")

iso_all_modes <- rbind(iso_5min_drive, iso_5min_walk)

otp_stop()
## [1] "SUCCESS: The process \"java.exe\" with PID 284 has been terminated."

Now I can draw a map of these isochrones on a map. I’ll use a background image from OpenStreetMap as a basemap (if you do this, always, always, always give attribution to the OpenStreetMap contributors). You’ll also want to add rosm.cache/ to your .gitignore file.

right_side <- st_bbox(iso_all_modes)$xmax
left_side  <- st_bbox(iso_all_modes)$xmin
top_side <- st_bbox(iso_all_modes)$ymax
bottom_side <- st_bbox(iso_all_modes)$ymin

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPL_libraries) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By car", "By foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")
## Loading required namespace: raster

OpenStreetMap has a few other basemap image options you can choose from, in addition to the default. You can see the list of options by typing rosm::osm.types() into your console

Here is cartolight:

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "cartolight", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPL_libraries) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By car", "By foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Here is stamenwatercolor:

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "stamenwatercolor", 
                      progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPL_libraries) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_discrete(name = "Area that is reachable within 5 minutes",
                      labels = c("By car", "By foot"),
                      type = c("gray", "black")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Here is stamenbw:

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "stamenbw", progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPL_libraries) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By car", "By foot")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Here is hotstyle:

ggplot(iso_all_modes) +
  annotation_map_tile(zoomin = 0, type = "hotstyle", 
                      progress = "none") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPL_libraries) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_discrete(name = "Area that is reachable within 5 minutes",
                      labels = c("By car", "By foot"),
                      type = c("gray", "black")) +
  theme_map() +
  labs(caption = "Basemap Copyright OpenStreetMap contributors")

Alternatively, you could just use the street network as a basemap.

ggplot(iso_all_modes) +
  geom_sf(data = cambridge_streets, color = "gray") +
  geom_sf(aes(fill = mode), alpha = 0.5) +
  geom_sf(data = CPL_libraries) +
  coord_sf(xlim = c(left_side, right_side), 
           ylim = c(bottom_side, top_side), expand = FALSE) +
  scale_fill_viridis_d(name = "Area that is reachable within 5 minutes",
                       labels = c("By car", "By foot")) +
  theme_map()

Calcuate and compare isochrone areas

I can use st_area() to calculate the area of each isochrone and visualize the relationship between the size of a walkshed and the size of a driveshed.

The function pivot_wider creates a separate column for each value of a specified variable (in this case, for each mode), so that each row represents a location (with three associated isochrones), rather than having each row represent an isochrone.

It’s not a flashy figure, but it is sort of interesting.

iso_areas <- iso_all_modes %>%
  mutate(area = st_area(iso_all_modes)) %>%
  st_set_geometry(NULL) %>%
  pivot_wider(names_from = mode, values_from = area) 

ggplot(iso_areas, 
       aes(x = as.numeric(walk), y = as.numeric(drive))) +
  geom_point() +
  scale_x_continuous(name = 
            "Area within a five-minute walking distance\nof a public library\n(square km)",
            breaks = breaks <- seq(10000, 130000, by = 20000),
            labels = breaks / 1000000) +
  scale_y_continuous(name = 
            "Area within a five-minute driving distance\nof a public library\n(square km)",
            breaks = breaks <- seq(0, 700000, by = 100000),
            labels = breaks / 1000000) +
  theme_bw()